Spin-orbital quantum liquid on the honeycomb lattice 
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In addition to low-energy spin fluctuations, which distinguish them from band insulators, Mott insulators 
often possess orbital degrees of freedom when crystal-field levels are partially filled. While in most situations 
spins and orbitals develop long-range order, the possibility for the ground state to be a quantum liquid opens 
new perspectives. In this paper, we provide clear evidence that the SU(4) symmetric Kugel-Khomskii model 
on the honeycomb lattice is a quantum spin-orbital liquid. The absence of any form of symmetry breaking - 
lattice or SU(N) - is supported by a combination of semiclassical and numerical approaches: flavor- wave theory, 
tensor network algorithm, and exact diagonalizations. In addition, all properties revealed by these methods are 
very accurately accounted for by a projected variational wave-function based on the :7r-flux state of fermions 
on the honeycomb lattice at 1 /4-filling. In that state, correlations are algebraic because of the presence of a 
Dirac point at the Fermi level, suggesting that the symmetric Kugel-Khomskii model on the honeycomb lattice 
is an algebraic quantum spin-orbital liquid. This model provides a good starting point to understand the recently 
discovered spin-orbital liquid behavior of Ba3CuSb209. The present results also suggest to choose optical 
lattices with honeycomb geometry in the search for quantum liquids in ultra-cold four-color fermionic atoms. 

PACS numbers: 67.85.-d, 71.10.Fd, 75.10.Jm, 02.70.-c 



I. INTRODUCTION 

The investigation of orbital physics in transition metal ox- 
ides has been recently boosted by the possibility to observe or- 
bital excitations with Resonant Inelastic X-ray Scattering[l]. 
This has been demonstrated in cases where the crystal-field 
splitting is strong enough to select a unique orbital configura- 
tion in the ground state and push the orbital excitations to high 
energy, hence to separate them from magnetic excitations. 
However, this is not the only possibility. If the electronic con- 
figuration of the transition metal ion is such that several orbital 
occupations are consistent with the crystal field environment, 
a situation referred to as orbital degeneracy [2], orbital fluctu- 
ations are likely to be much softer. In most cases known until 
recently, a cooperative Jahn-Teller distortion occurs, resulting 
in orbital order and gapped orbital excitations, but this needs 
not be the case a priori, and the search for situations in which 
orbitals remain fluctuating in the ground state has been very 
active over the past decade [3-9]. To which extent this oc- 
curs in the triangular system LiNi02 [10, 11] or in the spinel 
FeSc2S4 [12] is still debated [13]. Interestingly, a new candi- 
date has been recently put forward, Ba3CuSb209 [14], a Cu 
oxide that lives on a decorated honeycomb lattice in which no 
trace of orbital order could be detected. 

On the theory side, transition-metal oxides with orbital 
degeneracy are generally described by a Kugel-Khomskii 
model [15] in which spin and orbital degrees of freedom are 
coupled on each bond. A minimal model to investigate the 
possibility to stabilize an orbital liquid is the symmetric ver- 
sion of that model defined by the Hamiltonian 

^ = ^(2S, • S,- + ^)(2T,- . T,- + h (1) 

iij) 



where S/ is a spin- 1/2 operator and T/ a pseudo-spin 1/2 opera- 
tor that describes fluctuations of a two-fold degenerate orbital 
(a and b). Introducing the local basis |#) =|t | ) =\l a), 
|#) =1^^ b), I ) =\i b), the Hamiltonian exhibits the fufl SU(4) 
symmetry and can be written as = 2</j> Pij^ where Pij in- 
terchanges the states on sites / and j. The local basis states are 
often referred to as colors. In fact, the model is the straight- 
forward generalization of the SU(2) symmetric Heisenberg 
model for 5" = 1 /2 spins which, up to a constant and a fac- 
tor 2, has the same form when expressed with Pij operators. 

The first investigations of this model on various lattices 
have emphasized the role of 4-site plaquettes, the natural unit 
to build an SU(4) singlet [5, 16, 17]. The spontaneous for- 
mation of 4-site plaquettes has been proven for an SU(4) lad- 
der [17], and plaquette coverings have been argued to provide 
the relevant variational subspace for the ground state prop- 
erties of the SU(4) model on both the square and triangular 
lattices [5, 16], with possibly plaquette long-range order on 
the square lattice [18, 19]. In a variational study based on 
projected fermionic wave functions a gapless spin-orbital liq- 
uid state has been predicted in Ref. [8] on the square lattice. 
These previous conclusions have recently been challenged for 
the square lattice [20], for which spontaneous dimerization (as 
opposed to tetramerization) has been demonstrated on the ba- 
sis of state-of-the-art iPEPS (infinite projected entangled-pair 
state) simulations. In the dimerized phase, the dimers are an- 
tisymmetric states built out of 2 of the 4 colors, and they form 
a columnar state. Each dimer is preferentially surrounded by 
dimers built out of the other 2 colors, leading to long-range 
color order [20]. In the context of orbital degeneracy, this 
phase is likely to be ordered. Indeed, one of the states selected 
by the dimerization consists of alternating pairs of a and b or- 
bitals times a spin singlet. When coupled to the lattice, such a 



2 




FIG. 1. Summary of the main properties of the spin-orbital hquid of the SU(4) model on the honeycomb lattice, (a) Sketch of the local color- 
order. This pattern is the only one that respects the sequence of 4 colors along all the zigzag chains, whatever their orientation (horizontal, 
;r/3 or In 1 3). (b) Color structure factor of the Gutzwiller projected ;/r-flux state (VMC). The singular conical peaks are typical of algebraic 
correlations. 



state is expected to undergo a cooperative Jahn-Teller distor- 
tion that will stabilize orbitals a and b, hence to lead to orbital 
long-range order. 

In this paper, we consider the symmetric Kugel-Khomskii 
model on the honeycomb lattice. The first motivation is purely 
theoretical: since there are no 4-site plaquettes on this lattice, 
the ground state is unlikely to be a crystal of singlet plaquettes. 
The second one comes from experiments: the recent observa- 
tion of a spin-liquid behaviour in Ba3CuSb209 points to the 
honeycomb geometry as an outstanding candidate. 

The main result of the present investigation is summarized 
in Fig. 1: The SU(4) symmetric Kugel-Khomskii model is 
shown to be a quantum spin-orbital liquid with short-range 
color correlations that follow the pattern illustrated in the top 
panel, and strong evidence is provided in favor of an alge- 
braic spin-orbital liquid with typical conical singularities in 
the static structure factor, as shown in the bottom panel. 

To reach these conclusions, we have used a variety of an- 
alytical and numerical methods: Linear flavor- wave theory 
(LFWT), infinite projected entangled-pair states (iPEPS), ex- 
act diagonalization of finite clusters (ED), and a variational 
approach based on the Monte Carlo sampling of Gutzwiller 
projected fermionic wave-functions (VMC). Details about 
each method can be found in appendix A. These methods 
are complementary and shed lights on diff'erent aspects of the 
model: LFWT is a good starting point to test for lattice sym- 
metry breaking and color order. iPEPS is a variational ap- 
proach for infinite systems that has proven to be very suc- 
cessful to check the presence of any kind of long-range or- 
der [20, 21]. Exact diagonalizations reveal nearly exact in- 
formation on short length- scale properties and are extremely 
useful to benchmark other approaches. The variational Monte 
Carlo of fermionic wave-functions have proven to provide a 
remarkably accurate description of algebraic quantum liquids. 
As a first test, we compare in Fig. 2 the ground state energies 
of the various approaches. A number of conclusions can al- 



-0.75 
-0.8 
-0.85 

-0.95 



-1.05 



A 


VMC Majorana 7i-flux 




VMC 0-flux 


A 


VMC Majorana 0-flux 




-VMCTi-flux 




-iPEPS 2x2 unit cell 


-0- 


-IPEPS 4x4 unit cell 


□ 


ED 


T 


LFWT 



0.05 



0.1 0.15 
1/D or 1/N 



0.2 



0.25 



FIG. 2. Energy per site as a function of inverse bond dimension D 
(iPEPS) and as a function of inverse system size N (VMC and ED). 
We note that the LFWT energy is not variational. 



ready be drawn from this comparison. First of all, among all 
the Gutzwiller projected fermionic wave-functions we consid- 
ered, only one is really competitive, the wave-function based 
on the quarter-filled Fermi sea of standard fermions in the n- 
flux state (see Sec. I C for details). Its energy is much lower 
than that of the 0-flux state, as well as that of the half-filled 
Fermi sea of Majorana fermions with or tt flux, and these al- 
ternative fermionic wave-functions will not be considered any 
further^. Secondly, the agreement between iPEPS, ED and 
Gutzwiller projected n flux state is quite remarkable. This 
suggests that all these methods constitute appropriate descrip- 
tions within their range of validity.^ 



^ The variational energies (per site) of the tt-Aux, Majorana 0-flux, 0-flux, 
and Majorana 7r-flux states in the 96 site cluster are -0.895, -0.822, 
-0.763, and -0.755, respectively. 

^ LFWT leads to a very low energy. This is not a concern since the method 
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A. Absence of lattice symmetry breaking 

Lattice symmetry breaking, be it dimerization or plaquette 
formation, leads to bonds of different strengths. At the classi- 
cal level, which consists in minimizing the energy in the sub- 
space of product wave-functions of the form = Yli Wi)^ 
bonds are fully satisfied. Indeed, since the Hamiltonian of a 
bond Hij is a simple permutation, {i//ii//j\Hij\i//ii//j} = 
is minimal if neighboring states are orthogonal. On bipar- 
tite lattices such as the square or honeycomb lattices, it takes 
only 2 colors to achieve this, and the classical ground state for 
more than 2 colors is massively degenerate. This degeneracy 
however can be in principle lifted by zero-point fluctuations. 
The theory of harmonic fluctuations has been developed pre- 
viously, and it goes under the name of flavor-wave theory (see 
appendix A 1 for details). At the harmonic level, the energy 
of a bond takes the smallest possible value if the two colors of 
the bond are difl'erent from the other nearest neighbors. For 
the honeycomb lattice, this can be achieved for all bonds si- 
multaneously in an infinite number of ways. The configura- 
tion of Fig. 1(a) is the most symmetric one. In this configura- 
tion, all bonds have the same surrounding up to color permu- 
tations. Other configurations can be generated by exchanging 
the colors on a stripe of dimers (see Fig. 3(a-b)), leading to a 
degeneracy of order 2 ^ since, once a direction has been cho- 
sen, this can be done independently on all dimer stripes. In 
all configurations, the energy is the same on all bonds. This 
is a first hint that, by contrast to the square lattice, the lattice 
symmetry is not broken for the honeycomb lattice. 

The same family of degenerate ground states is obtained 
with iPEPS if a unit cell with 4x4 = 16 diff'erent tensors 
and a small bond dimension Z) = 2 are used. Upon increas- 
ing D, more quantum fluctuations are taken into account, and 
the symmetric state of Fig. 3(a) is stabilized. In this state, all 
bonds have the same energy (see Fig. 3(c)). To test how robust 
this conclusion is, we have challenged it by performing iPEPS 
simulations using a 2 x 2 unit cell with only four difl'erent ten- 
sors. This leads to a dimerized state with two types of dimers 
A and B, which can be distinguished by their dominant colors, 
and difl'erent inter- and intra-dimer bond energies (Fig. 3(d)). 
However, unlike on the square lattice, this dimerization van- 
ishes in the infinite D limit, as shown in Fig. 4(a) where the 
diff'erence in bond energies, AEi, = max(£/,)-min(£'/^), is plot- 
ted. Thus, both low-energy states found with iPEPS preserve 
the lattice symmetry in the large D limit. 

We also tested with VMC the stability of the 7r-flux state 
towards the dimerization instability shown in Fig. 3(d) by 
strengthening/weakening the hopping amplitude of the bonds, 
with the conclusion that the variational energy is minimal in 
the absence of any dimerization. Similarly, we found no in- 
dication toward quadrumerization (SU(4) singlet formation), 
where the hoppings connected to say red sites (forming a 'tri- 
pod') in Fig. 1(a) are modified. 




FIG. 3. Examples of states obtained with LFWT (a-b) and iPEPS for 
small bond dimension Z) (c-d). (a) Most symmetric configuration, (b) 
Configuration obtained from the most symmetric one by exchanging 
colors on a stripe, (c) Color-ordered state with one dominant color 
per site, obtained with a 4 x 4 unit cell and a bond dimension D - 6. 
(d) Dimer-Neel ordered state obtained with a 2 x 2 unit cell (shaded 
in blue) and D - 6. Both the color order and the dimerization vanish 
in the infinite D limit. The pie charts show the local color density on 
each site and the thickness of a bond is proportional to the square of 
the energy on the corresponding bond. 



Finally, in Fig. 4(b) we show the ED results for the con- 
nected bond energy correlations, which is a way to detect 
dimerization and other bond energy pattern formation tenden- 
cies. The correlations decay quite rapidly with distance, mak- 
ing dimerization or other patterns unlikely. 

Thus, all methods consistently point towards a state which 
does not break the lattice symmetry. 



B. Absence of SU(4) symmetry breaking 

The color-ordered states predicted by LFWT and iPEPS 
with a small bond dimension (Fig. 3) break the SU(4) sym- 
metry. Here we show that higher-order quantum fluctuations 
destroy this color order, i.e. that in the ground state the SU(4) 
symmetry is in fact unbroken. 

In Fig. 4(c) we present the iPEPS result for the local ordered 
moment. 



^ 



(2) 



is not variational, but this is an indication that, in the present case, higher 
orders are likely to be important. 



where S^, = \a}{j3\ are the generators of SU(4) and a, 13 run 
over the four difl'erent flavors. A finite m implies that the 
SU(4) symmetry is spontaneously broken in the thermody- 
namic limit. For both low-energy states found with iPEPS the 
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FIG. 4. Various correlation functions obtained with iPEPS, ED 
and VMC. (a) Difference in bond energies obtained with iPEPS in 
the state shown in Fig. 3(d). It is strongly suppressed with increasing 
bond dimension D and vanishes in the infinite D limit, (b) Connected 
bond energy correlations (PijPu) - (Pij)(Pki) calculated with ED in 
the ground state of the N = 24 sample. The black bond denotes 
the reference bond. Solid blue (dashed red) bonds stand for positive 
(negative) correlation functions, and the width is proportional to the 
absolute value of the correlation function, (c) Local ordered moment 
m obtained with iPEPS as a function of inverse bond dimension. It 
vanishes in the infinite D limit for both low-energy states shown in 
Fig. 3(c-d). (d) Spin correlation function in real space, as calculated 
from ED (right) and VMC (left) for a 24- site cluster. The area is 
proportional to (Pq/) - 1/4, where is the index of the central site, 
and / labels the sites in the 24- site cluster. The color keeps track of 
the sign (blue for positive, red for negative). 
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FIG. 5. Properties of the ;7r-flux state state, (a) Sketch of the gauge 
used to implement the ;7r-flux state: the hopping amplitude is positive 
on solid blue bonds, negative on dashed red bonds. The primitive 
unit cell (dark magenta) contains four sites, the hexagonal unit cell 
eight sites, (b) Brillouin zones and high symmetry points. The red 
circles indicate the position of the Dirac-nodes slI sd = + V3? to 
which the Fermi surface reduces at quarter filling in the ;r-flux state. 
The orange, outermost hexagon shows the extended Brillouin zone 
of the triangular lattice (including sites at the centers of the hexagons 
in the honeycomb lattice), the structure factor is maximal and has a 
cusp at Ma = (n,7T/ V3) and the symmetry related points, is given 
by (4;r/3, 0), K is (2;r/3, 2;r/3 V3). (c) The two-fold degenerate band 
structure of the ;7r-flux state in the reduced Brillouin zone of 8-site 
hexagonal unit cell. 



liquid. 



C. Algebraic spin-orbital liquid 



local ordered moment is strongly suppressed with increasing 
bond dimension, and most likely vanishes in the large D limit. 
This shows that quantum fluctuations, which are systemati- 
cally taken into account by increasing D, eventually destroy 
the color order so that the SU(4) symmetry is restored. This is 
in contrast to the same model on the square lattice [20] where 
the local ordered moment was found to remain finite in the 
infinite D limit. 

Consistent results for the flavor correlation function are 
obtained with ED and VMC for the 7r-flux state shown in 
Fig. 4(d), which is decaying rapidly with increasing distance, 
indicating absence of long range order. The very good qualita- 
tive and quantitative agreement between the ED and the VMC 
results provides substantial evidence that the 7r-flux state cor- 
rectly describes the short range physics of the ground state of 
the Hamiltonian (1). In the next section we will show that 
the decay predicted by VMC is algebraic, i.e. that the state 
described by this wave function is an algebraic spin-orbital 



A standard way to describe spin liquids for SU(2) models 
is based on the fermionic representation of the spin operators 
[22-26] using a variational wave function [27, 28], where the 
multiply occupied sites are projected out from a suitable cho- 
sen, noninteracting Fermi-sea. In the generic case, the Fermi- 
sea has a finite Fermi surface. However this is not the only 
possibility. For the SU(2) Heisenberg model on the square 
lattice, Marston and Afiileck have shown that introducing a tt- 
flux per elementary plaquette leads to the formation of Dirac- 
nodes [29]. At half filling, the Fermi surface of this 7r-flux 
state shrinks to points, and its energy is lower than that of 
the state with equal hopping amplitudes and a finite Fermi 
surface. In such a spin-liquid the structure factor is singular 
at momenta related to the diflTerence between Fermi-points, 
leading to the algebraic decay of spin correlations. In one 
dimension, this type of approach leads to an accurate descrip- 
tion of the algebraic decay of the correlations for the SU(2) 
case [30], and as well as for SU(4) using the representation 
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On the honeycomb lattice, a Dirac-node is already present 
at the middle of the band without any flux, and the Fermi sur- 
face reduces to points at half filling. So the zero flux state 
would be a good starting point to describe an algebraic spin 
liquid for the SU(2) Heisenberg model [31-33]. However, for 
the SU(4) Heisenberg model, the band must be quarter-filled, 
and the equivalent of the Affleck-Marston approach requires 
to have a Dirac node at the Fermi energy of the quarter-filled 
system. It turns out that this is achieved in the 7r-flux state, as 
shown in Fig. 5(c). As for the square lattice, this state leads 
to a lower energy than the zero-flux state, as already stated 
above. 

Starting from the noninteracting wave function, with a band 
populated up to the Dirac-node at Sd = -^[?)t for any of 
the four flavored fermions, we implemented the Gutz wilier- 
projection using a variational Monte-Carlo (VMC) sampling. 
The energy of this wave function, E = -0.894 per site, com- 
pares remarkably well with that of iPEPS (see Fig. 2), espe- 
cially considering that no variational parameter was used. Let 
us also mention that the state (and the ones related by symme- 
try) shown in Fig. 1(a) has the maximal weight in the varia- 
tional wave function. 

To investigate the physics of this wave function, we have 
calculated the spin-spin correlation function as a function of 
distance. The results clearly demonstrate an algebraic decay 
{Pij - 1/4) ~ Ir/ji"^, with an exponent a between 3 and 4, 
as shown in Fig. 6. If one considers the honeycomb lattice 
as built from zigzag chains, these correlations correspond to 
even distances along one of the zigzag chains, and the expo- 
nent should be compared to that of the dominant correlations 
with wave vector 7r/2 of a single chain [34]. This exponent 
is equal to 3/2, a number actually very accurately reproduced 
by VMC. So color-color correlations decay faster on the hon- 
eycomb lattice than on a chain, but still algebraically. This is 
a rather peculiar situation in view of the standard paradigms: 
the development of long-range order, as in weakly coupled 
SU(2) chains in square geometry, or the spontaneous forma- 
tion of local singlets and exponentially decaying color-color 
correlations, as e.g. in the SU(4) ladder [17]. 

This Gutz wilier projected tt flux state is actually a proto- 
typical wave function for a phase which should be called al- 
gebraic spin-orbital liquid in the present context, in analogy 
to the algebraic spin liquids which have been discussed in 
the spin liquid literature [22, 24, 25, 35]. These states are 
characterized by the algebraic decay of a number of correla- 
tions, with wave vectors corresponding to difl'erences between 
the Dirac cone locii. In our case for example the algebraic 
spin-orbital correlations are modulated with wave vector M in 
the standard honeycomb Brillouin zone, and this corresponds 
to the distance between two Dirac cones in the tt flux state. 
Another important aspect of this wave function is that it has 
been shown that it can describe an extended phase in param- 
eter space and not just an unstable fine-tuned point [24, 36]. 
Upon adding perturbations of suitable strength, many diff'erent 
phases can be found in the vicinity of an algebraic spin-orbital 
liquid [24], making the present model an interesting starting 
point for further explorations of exotic phases in spin-orbital 
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FIG. 6. The algebraic decay of the correlation along a zig-zag chain 
in the honeycomb lattice for the Gutz wilier projected ;r-flux state in 
different clusters (every second site is shown, 6 is the distance along 
the chain). In the inset we compare the correlations of the same tt- 
flux state with the correlations of a 300 site long one-dimensional 
chain (we projected the one-dimensional quarter-filled Fermi- sea). 
While the periodicity of four is visible in both cases, the correlations 
decay much faster in the two-dimensional honeycomb lattice. 



systems. 



II. CONCLUSIONS 

In conclusion, the results reported in this paper provide very 
strong evidence that the SU(4) symmetric Kugel-Khomskii 
model is a quantum spin-orbital liquid, and build a case in 
favor of an algebraic spin-orbital liquid. Clearly, the present 
results do not allow to exclude that the quantum spin-orbital 
liquid is of another type, for instance some kind of resonating 
valence-bond (RVB) liquid with resonances between 4-site 
cluster singlets, but our results suggest that the correspond- 
ing gap would be quite small. In particular, experience with 
highly frustrated magnets have shown that, as good as its en- 
ergy might be, a fermionic variational wave function might 
fail to capture the correct low-energy physics. This seems 
for instance to be the case for the SU(2) Heisenberg model 
on the kagome lattice, for which the variational energy of the 
algebraic spin liquid wave function [25] is close to the best 
numerical estimates [37, 38], yet DMRG results have given 
strong evidence in favor of a gapped Z2 spin liquid [37]. In 
the present case, the projected fermionic wave function has 
not just been tested for its energy, but correlations have been 
shown to be in remarkable agreement with ED up to interme- 
diate distances. So we believe that the case is strong, but not 
closed. 

In any case, the fact that the ground state is a quantum spin- 
orbital liquid is quite firmly established. This result is quite 
interesting in view of the liquid behavior reported recently 
in Ba3CuSb209. A detailed comparison will clearly require 
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some adaptation of the present model, in particular to take 
into account the additional magnetic Cu sites present in the 
system on top of the honeycomb lattice. Still, the absence of 
orbital order reported in the present paper is consistent with 
the experimental results. 

Finally, we note that the SU(N) Heisenberg model is a 
rather accurate effective model for the 1/A/^-filled Mott insu- 
lating phase of alkaline-earth metal atoms with N internal de- 
grees of freedom loaded in an optical lattice [39-41]. Cur- 
rently the main issue in that field is to reach low enough en- 
tropies to observe correlations typical of long-range order, but 
the next step will definitely be to realize exotic quantum states. 
In that respect, the = 4 case on the honeycomb lattice ap- 
pears to be a very strong candidate. 
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Appendix A: Methods 

1. Linear flavor- wave theory 

The linear flavor- wave theory (LFWT) is a method to treat 
harmonic quantum fluctuations on top of a mean-field (or 
Hartree) solution based on a site-factorized variational wave 
function [42, 43]. It starts from a Schwinger boson repre- 
sentation of the SU(4) operators with 4 types of bosons. In 
a mean-field ground state, each site has a well defined color. 
The corresponding boson is assumed to condense, and the re- 
sulting Hamiltonian is a bosonic quadratic form. On a nearest- 
neighbor bond with color a on site / and color JS on site j, it is 

given by "Hfw = A].A. . - 1 , with A]. = bl .-\- br, ■ where b^ and b 

u u u p,i 

are bosonic operators. In a given mean-field ground state, the 
LFWT Hamiltonian is the sum of independent Hamiltonians 
that describe the motion of bosons on the connected clusters 
spanned by pairs of colors. The zero-point energy per bond 
tends to increase with the cluster size, and is minimal on a 
2- site cluster, when the ground state energy of the Hamilto- 
nian is equal to -1, i.e., there is no zero-point contribution to 
the energy. By contrast, larger clusters lead to finite frequen- 
cies, hence to strictly positive contributions to the zero-point 
energy. 

2. Infinite projected entangled-pair states (iPEPS) 

A projected entangled-pair state (PEPS) [44, 45], also 
called tensor product state, is a variational ansatz where the 
wave function of a two-dimensional system is efficiently rep- 
resented by a product of tensors, with one tensor per lattice 
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FIG. 7. (a) The honeycomb lattice is mapped to a square lattice with 
brick-wall structure. There is no Hamiltonian term between sites 
connected by a dotted line. A square lattice iPEPS is used for this 
lattice, where we choose the bond dimension along the dotted lines 
as Z) = 1 . (b) Ordered moment as a function of inverse x which 
controls the accuracy of the contraction of the iPEPS. The values of 
m for different bond dimensions D only slightly depend on x- 



site. It can be seen as a two-dimensional generalization of 
matrix product states - the class of variational states underly- 
ing the famous density matrix renormalization group (DMRG) 
method [46] . On the square lattice each tensor T^.^^ has a phys- 
ical index p carrying the local Hilbert space of a lattice site 
with dimension J, and four auxiliary bond indices /, j, k, I with 
dimension D which connect to the four neighboring tensors. 
Thus, each tensor consists of JZ)^ variational parameters, and 
by changing the bond dimension D the accuracy of the ansatz 
can be systematically controlled. A Z) = 1 PEPS simply corre- 
sponds to a site factorized wave function (a product state), and 
upon increasing D quantum fluctuations can be systematically 
taken into account. 

An infinite PEPS (iPEPS) [47, 48] consists of a unit cell 
of Lx X Ly = Nt tensors which is periodically repeated in 
the lattice to represent a wave function in the thermodynamic 
limit. We use the iPEPS method developed for the square lat- 
tice described in Refs. [20, 47, 49] to simulate the model on 
the honeycomb lattice by mapping it onto a brick- wall lattice 
as illustrated in Fig. 7(a). The bond dimension of the auxiliary 
bonds connecting two sites which do not directly interact (dot- 
ted lines) can be chosen sls D = 1 . This ansatz is equivalent to 
an iPEPS with only three auxiliary indices on the honeycomb 
lattice. The advantage of this mapping is that the codes devel- 
oped for the square lattice only require minor modifications to 
simulate models on the honeycomb lattice. 

Describing the iPEPS method in full detail is beyond the 
scope of this paper and we therefore only mention the most 
important technicalities with corresponding references and 
details on the simulation parameters in the following. For an 
introduction to PEPS and iPEPS we refer to Refs. [45, 49]. 

The optimization of the tensors (i.e. finding the best vari- 
ational parameters) is done through imaginary time evolu- 
tion, first with the so-called simple update (see Refs. [49, 50]) 
which is equivalent to the time-evolving block decimation 
method in one dimension [51, 52]. The solution is used as 
an initial state for an imaginary time evolution using the full 
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update [47, 49], which is computationally more expensive, 
but more accurate than the simple update, since the full wave 
function is taken into account for the truncation of a bond 
index. We use a second order Trotter- Suzuki decomposition 
with time steps down to dr - 10"^. For large values of Z) a 
larger time step dr = 10"^ is used, where the estimated Trotter 
error is small compared to symbol sizes (e.g. below 0.5% for 
the ordered moment m (2)). 

Expectation values of observables can be computed by con- 
tracting the tensor network, i.e. by computing the trace of 
the product of all tensors. For the approximate contraction of 
the iPEPS we use the corner transfer matrix method described 
in Refs. [49, 53, 54]. The accuracy of this contraction can 
be controlled by the so-called boundary dimension where 
we used values up to ;^ = 500 (typically up to ;^ = 350) for 
large D. Observables like the energy and the ordered moment 
are extrapolated in with an extrapolation error being small 
compared to symbol sizes. An example is shown in Fig. 7(b) 
for the ordered moment m (2). 

To increase the efficiency of the method we implemented 
a Zq symmetry in the tensors (see Refs. [55, 56]) which is a 
discrete subgroup of SU(4). This implies that the tensors have 
a block structure (similar to a block diagonal matrix) which 
reduces the computational cost considerably. 

Since iPEPS is an ansatz for an infinite system, symme- 
tries such as SU(4) or translational symmetry may be sponta- 
neously broken. In order to test diff'erent types of translational 
symmetry breaking we compared the variational energies ob- 
tained with diff'erent unit cell sizes. We found two competing 
low-energy states with unit cell sizes 2x2 and 4x4, shown in 
Fig. 3, which have similar energies for large bond dimension. 
We note that broken symmetries can be an artifact of a finite 
bond dimension D, and that the symmetry can be restored if D 
is sufficiently large. In other words, a classical or a low entan- 
glement solution (small D) might exhibit order, but this order 
can be destroyed by quantum fluctuations which are systemat- 
ically included upon increasing D. Therefore, it is important 
to study order parameters as a function of D to verify if they 
are finite in the large D limit. 



3. Exact diagonalization 

We have performed exact diagonalizations of the Hamil- 
tonian (1) for selected finite size samples of up to A/^ = 24 
sites. The energies of the samples with 8,14,18 and 24 sites 
are reported with star symbols in Fig. 2. Note that only the 
samples with 8 and 24 sites can form a SU(4) singlet ground 
state, and this explains the relatively high energy per site for 
the samples with 14 and 18 sites. Note also that due to com- 
putational limitations we were only able to calculate the en- 
ergy and eigenfunction of the ground state in the completely 
symmetric representation of the spatial symmetry group of the 
N = 24 sample (the Hilbert space in this symmetry sector con- 
tains 4' 008 '4 17 '658 states, including a cyclic color permuta- 
tion symmetry). Since this sector is the absolute ground state 
for the A/^ = 8 sample, we expect this sector to host the ground 
state for A/^ = 24 as well. 



4. Fermionic Variational Monte Carlo 

The variational wave function for the algebraic spin-orbital 
liquid is a Gutz wilier projected noninteracting Fermi sea at 
quarter filling: 



(Al) 



where denotes the position of the /-th fermion with color a, 
and the summation is over the all N\/((N/ 4)1)^ possible dis- 
tributions of the fermions so that each site is single occupied 
(i.e. {7^} n {f} is an empty set for a 9^ JS). The weight of each 
configuration is given by the 



(A2) 



Slater-determinant, where ^k(j) is the amplitude of the 
fermion at site j in the ^-th eigenfunction of the hopping 
Hamiltonian 



(ij) a=l 



(A3) 



The expectation values of operators with this variational wave 
function are evaluated by a Monte Carlo sampling [27, 28]. 
The variational wave function with Majorana fermions is de- 
scribed in detail in Ref. [8]. 

The diff'erent trial states correspond to diff'erent choices of 
thexij hopping parameters. In the ;r-flux state \xij\ = 1 and 
the phases are chosen so that an electron hopping around the 
hexagon picks up a minus sign, FlLi^^.^+i ~ where the 
product is over the bonds of a hexagon. This can be realized in 
many ways. Here we choose real hopping amplitudes, where 
every hexagon has a single bond with a Xij = - 1 . while the 
rest of the bonds have -hi, as shown in Fig. 5(a). Furthermore, 
we allow for antiperiodic boundary conditions when a degen- 
eracy of quarter filled Fermi sea needs to be removed for a 
given cluster. 

We considered two families of finite size clusters with the 
full symmetry of the honeycomb lattice: (i) clusters with 
N = lilnf sites defined by the lattice vectors gi = ( V3, 0)^ 
and g2 = ( V3/2, 3/2)^, where n is an integer (A^ = 72, 200, 
392, and 648 site clusters); (ii) clusters with N = 6(2nf sites 
with gi - (3 V3/2, 3/2)^ and g2 = (0, 3)^, (A^ =24, 96, 216, 
384, 600). 

We use Monte Carlo sampling of the projected wave func- 
tion 1^) to evaluate physical quantities. The elementary step 
is the exchange of two randomly chosen fermions with dif- 
ferent colors. To speed up the convergence we used impor- 
tance sampling, with acceptance ratios defined according to 
the Metropolis algorithm: the new configuration is always ac- 
cepted if its weight Wnew is higher than the weight Wom of the 
old configuration, but for Wnew < >^oid it is accepted only with 
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a probability Wnew/>^oid- The configurations are thus repre- 
sented with a probability p{j} proportional to their weight in 
the wave function: 



P{j} ^ 



(A4) 



where W{j} = Yla^i ^{f} denotes the coefficient of \{j}) = 
'^t^i^JiJi " ' ill th^ projected fermionic state |^), see 
Eq. (A2). We set the number of elementary steps between two 
measurements large enough to avoid autocorrelation eff'ects. 
See Table I for details. We also performed a binning analysis 
as a further test for the independence of the measurements. 
The statistical error for diff'erent bin sizes did not show any 
change, verifying once more that the sampling distances were 
large enough. Furthermore, for each system we made several 
(5-10) runs with randomly chosen starting configurations to 
independently verify the error bars obtained from the binning 
analysis. 



N 


^a.c. 


An 


ratio 


number of measurements 


24 


22 


1000 


45.5 


10^ 


72 


150 


1000 


6.7 


10^ 


96 


260 


2000 


7.7 


10^ 


200 


970 


5000 


5.1 


2- 10^ 


216 


1080 


5000 


4.6 


2- 10^ 


384 


2920 


20000 


6.8 


2- 10^ 


392 


3340 


20000 


6.0 


2-10^ 


600 


7080 


40000 


5.6 


10^ 


648 


8100 


40000 


4.9 


10^ 



TABLE I. Ta.c. autocorrelation times for the two site correlation func- 
tions compared to the number of elementary step between two mea- 
surements (An). 

We measured diagonal and off'-diagonal operators. The 
spin-spin correlation function, the average of the off'-diagonal 
Pkj can be expressed using the diagonal operator (where 
^ is the occupation number on site k for the fermion of color 

as 



supposing that the ground state is a singlet wave function - as 
it is the case when the hopping Hamiltonian is independent of 
the colors. The measurement of the diagonal (^^) correla- 
tion functions is quite simple using the importance sampling. 



MC 



N({kJ}c{f}) 
Nmc 



(A6) 



where N({k, 1} c {f}) denotes the number of times both k and / 
sites were occupied with JS fermion among the A^mc measured 
configurations. 

With a little more eff'ort one can directly calculate the 
off'-diagonal {Pk,i) correlation functions as well. Using the 
fermionic representation for the Pk,i exchange operator, the 
convenient form that follows the convention of the fermion 
ordering in the wave function is given as 

pk,i = Yj ^a(k)s-^(i) = - ^ flm^(% ma (0 (A7) 

aj3 aj3 

In this case one follows the same importance sampling as be- 
fore, although the measurement itself is more complicated: 



^uuj} '^oi'^iii ((ilin,/ 1{ J}) 



MC ^{j} 
{J\mc 



(AS) 



(A5) 



where {/} is the configuration that leads to {7} by exchanging 
the color of fermions on sites k and /, and the sum in the last 
equation is over the measured configuration [ jjuc- Following 
Eq. (A7), the sign Skj({ j}) is 1 if the colors of fermions on 
sites k and / in the configuration {j} are the same, and -1 if 
the colors are diff'erent. We have explicitly verified that the 
Eq. (A5) holds. 

Similarly one can calculate (PijPki) as well, here for each 
{7} configuration \{j}) = PijPkiWf})- For the sign one should 
take Skj({ j})sij({ j}), here we assumed that the sites ij,k, and / 
are all diff'erent. 
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